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Abstract 

In a recent paper [T] the scattering and transport of excess electrons in liquid argon in the hydro- 
dynamic regime was investigated, generalizing the seminal works of Lekner and Cohen with 
modern scattering theory techniques and kinetic theory. In this paper, the discussion is extended 
to the non-hydrodynamic regime through the development of a full multi-term space-time solution 
of Boltzmann’s equation for electron transport in gases and liquids using a novel operator-splitting 
method. A Green’s function formalism is considered that enables flexible adaptation to various 
experimental systems. The spatio-temporal evolution of electrons in liquids in the hydrodynamic 
regime is studied for a benchmark model Percus-Yevick liquid as well as for liquid argon. The 
temporal evolution of Franck-Hertz oscillations are observed for liquids, with striking differences 
in the spatio-temporal development of the velocity distribution function components between the 
uncorrelated gas and true liquid approximations in argon. Transport properties calculated from the 
non-hydrodynamic theory in the long time limit, and under steady-state Townsend conditions, are 
benchmarked against hydrodynamic transport coefficients. 
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I. INTRODUCTION 


The study of electron transport in dense gases and liquids is of interest in nnderstanding 
the fundamental microscopic scattering processes involved, and to technological applications 
including liquid-state electronics [1], high-energy particle detectors [HHi, plasma medicine 
[SHE] and radiation dosimetry [EHIS]. For these technologies to reach their full potential 
reqnires a detailed understanding of the fnll spatio-temporal behavionr of electrons in dense 
gases, liquids, and other bio-structures, typically under non-equilibrinm conditions. 

In a recent paper [T] the transport of excess electrons in liquid argon was considered from 
ab initio liqnid phase cross-sections calcnlated nsing the Dirac-Fock scattering eqnations. 
The approach detailed in the seminal works by Lekner and Cohen Els] has been revisited 
with modern scattering theory techniques, where the original treatment was extended to 
consider multipole polarizabilities and a non-local treatment of exchange. With an increase 
in density, several important density effects become significant, most notably (i) the coherent 
scattering from multiple scattering centres, (ii) the screening of the long range polarization 
potential dne to induced multipoles in the bnlk, and (iii) the contribution of the bulk to the 
effective potential experienced by the electron. Transport coefficients snch as drift velocities 
and characteristic energies calculated in the hydrodynamic regime with our hydrodynamic 
multi-term Boltzmann equation solution were in good agreement with swarm experiment 
measurements in both gas- and liquid-phase argon na. In this work we extend the discussion 
to an investigation of liqnid state in the non-hydrodynamic regime, using the same electron- 
argon potentials and cross-sections presented in [T]. 

The solution of the the fnll temporal-, spatial- and energy-dependent Boltzmann eqnation 
is formidable, both mathematically and compntationally. Historically, the majority of kinetic 
theory investigations have focnsed on the hydrodynamic regime where spatial gradients are 
small, and have considered increasingly complex space- and time-dependent hydrodynamic 
behaviours and field configurations (see reviews [IBH22])- In situations where the hydrody¬ 
namic regime is not applicable, the space-time dependence of the phase space distribution 
function cannot be projected onto the number density and a density gradient expansion 
is no longer valid. Instead the confignration-space dependence of the Boltzmann eqnation 
mnst be treated on eqnal footing with the energy-space dependence, which makes for a dif¬ 
ficult problem even for simple geometries |2SH23- It is no surprise that systematic studies 
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of non-hydrodynamic phenomena lag behind their hydrodynamic counterparts. The proto¬ 
typical example of non-hydrodynamic phenomena is the Frank-Hertz experiment [SSIEZI, 
which helped lay the foundations for quantum and atomic physics. Extensive theoretical 
studies of non-hydrodynamic electron phenomena have been performed including field free 
spatial relaxation [28] , and spatial relaxation in the presence of uniform [2D]-[32| , non-uniform 
[33] and periodic electric fields |Ml - l36| . Similar kinetic studies on the spatial relaxation of 
electrons in uniform and spatially periodic fields have been performed by Golubovskii et 
al. [HTHIH] . Li and co-workers have considered arbitrary electric and magnetic field con¬ 
figurations with a multi-term analysis naiiiiiii. Solution of the Boltzmann equation for 
electrons including both the space and time dependence have also recently been performed 
|i3l - H5] . however these authors restricted their calculations to a two-term approximation in 
Legendre polynomials in order to make the problem computationally feasible. Limitations 
of the two-term approximation for molecular gases are well known m- Prior to [T|, all 
studies of electron transport in liquids were in the hydrodynamic or spatially homogeneous 
regimes, and restricted to the two-term approximation. 


In this study, we present a full multi-term space-time dependent solution of Boltzmann’s 
equation, capable of handling highly non-equilibrium electron transport in dilute gases, dense 
gases and liquids under non-hydrodynamic conditions. To our knowledge, this is the first 
time such a complete solution of Boltzmann’s equation has been developed. In addition, by 
solving for the spatio-temporal evolution of the Boltzmann equation Green’s function, the 
technique is quite general in its application, enabling various experimental configurations 
(temporal and spatial initial and boundary conditions) and practical devices to be modelled 
from a single solution. This work extends the Boltzmann equation framework to applications 
and accuracies comparable to those achieved using the Monte-Garlo simulations of Petrovic, 
Dujko and co-workers [171 - 11^ . 


We begin the paper in Section [IT] with a brief overview of the multi-term solution of 
Boltzmann’s equation for electrons in structured materials, such as liquids. We then detail 


our operator splitting treatment of the space, time and energy dependence in Section |TTT| In 
Section we present solutions for a model hard-sphere liquid system with a Percus-Yevick 
structure factor used to simulate a prototypical liquid with realistic pair correlations. A 
simple inelastic channel is included to induce periodic oscillatory structures (an idealized 
verion of the well known Frank-Hertz experiment 123) which can act as a non-hydrodynamic 
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benchmark. Lastly, in Section |V] we investigate the temporal and spatial evolution of the 
phase-space distribution for electrons in liquid argon, using microscopic cross-sections which 
have been derived previously |Tj . The issues with treating liquid systems as gaseous systems 
with increased density, and the implications for various applications including liquid argon 
time projection chambers, are highlighted. 


II. THEORY 

The fundamental kinetic equation describing the evolution of an electron swarm in a 
gaseous or liquid background medium subject to an electric field, E, is the Boltzmann 
equation for the phase-space distribution function / = / (r, v,t) [HO] : 


d d eEd 

dv~^ m &v ) 




-JU), 


( 1 ) 


where t is time, and r, v, e and m are the position, velocity, charge and mass of the electron 
respectively. The collision operator J{f) accounts for interactions between the electrons and 
the background material, and describes the effect of collisions on the distribution function 
at a fixed position and velocity. In essence the Boltzmann equation represents an extension 
of the continuity equation to phase-space. A solution of equation Q for the distribution 
function yields all relevant information about the system. Macroscopic transport properties 
including mean energy and drift velocity can then be found via averages over the ensemble, 
as detailed in Section ITU FI 

The starting point for most modern solutions of the Boltzmann equation is the decom¬ 
position of the angular part of the velocity dependence of equation 0 in terms of spherical 
harmonics If there is a single preferred direction in the system, e.g. due to an electric 
field in plane parallel geometry, then the angular dependence of the velocity component 
can be adequately described by a simpler expansion in terms of Legendre polynomials. 
For the plane-parallel geometry considered in this work, /(v,r,t) —?• f {v, z, where 
/i = V • E = cos X, such that 


/(v,r,t) = ^fi {v,z,t)Pi{fi), 
1=0 


( 2 ) 
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where Pi is the /-th Legendre polynomial. Upon substituting the expansion Q into equa¬ 
tion 0 and equating the coefficients of the Legendre polynomials results in the following 
system of coupled partial differential equations for the /p 


dh 

dt 


p=zbl 



u''^-^ + eE 

oz 


1 Q U + 1 


fl+p 


-Ji(fi), (3) 


for I = 0,1, 2,..., oo, where U = Ji is the Legendre decomposition of the collision 

operator, and 


. (+ 1 ) ^ I . (- 1 ) ^ {I + 1 ) 

^ (2/-1)’ ^ (2/+ 3)' 


(4) 


Equation Q represents an infinite set of coupled partial differential equations for the expan¬ 
sion coefficients, //. In practice, one must truncate the series Q at a sufficiently high index, 
I = ^max- The history of charged particle transport in gases and liquids has been dominated 
by the ‘two-term approximation’ , i.e., where only the first two terms have been included. 
The assumption of quasi-isotropy necessary for the two-term approximation is violated in 
many situations, particularly when inelastic collisions are included lai or when higher order 
moments are probed [Tj. Such an assumption is not necessary in our multi-term formalism. 
Instead, the truncation parameter is treated as a free parameter that is incremented 
until some convergence criterion is met on the distribution function or its velocity moments. 

In order to solve equation 0 we require the collision operators for all of the relevant 
collisional processes, and their representations in terms of Legendre polynomials, Ji. If we 
assume that the neutral background material is at rest and in thermal equilibrium at a 
temperature Tq, then the collision operator is linear in the swarm approximation. Below we 
detail the specific forms of the collision operator for particle-conserving elastic and inelastic 
collisions employed. A further expansion of each collision integral with respect to the ratio 
of swarm particle mass to neutral particle mass, m/mo, has been performed. Because this 
ratio is small for electrons in argon, only the leading term of this expansion was taken into 
account for each collision process and in each equation of the system Q . The total collision 
operator can be separated for each of the different types of processes, e.g. 


Ji 


j; 


el 


I rexc 
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where Jf and are the elastic and inelastic collision operators, respectively. For dense 
mediums and low electron energies, the de Broglie wavelength of the electron is comparable 
to the average inter-particle spacing ^ tTq . The electron wave is then scattered coherently 
from multiple scattering centres in the medium. Cohen and Lekner [3] showed how to account 
for coherent scattering using a two-term approximation, which has since been extended to a 
multi-term regime [501 [52]. The Legendre projections of the elastic collision operator in the 
small mass ratio limit were shown to be: 


Jf ifi) = 


‘2mjT-J d 
rriQ^ dU 

vfiumu) 


UivfiU) (A + faToH) 


I = 0 , 

I > 1 . 


( 6 ) 


Here, uiq is the mass of the background neutral, kb is Boltzmann’s constant, and 


ufiU) =no{^) 
' m ) 


2vr / a(t/,x) [1 - Pi{cosx)] sinydy ) , 


(7) 


are the usual binary collision frequencies. A collision frequency, p, is related to the cor¬ 
responding cross section, cr, via p = Uq cr({7). The pf{v) in equation (6) are the 

structure-modified counterparts to pf{U), i.e.. 


i)f{U) = no{—y (27: 


S(Vx) [1 - a (cosy)] sinydy ) , 
where E(f7, y) is the effective differential cross-section 


m / \ JO 


( 8 ) 


S(17,y) = a{U,x) S Qv^^sin^ , (9) 

which then accounts for coherent scattering effects through the static structure factor, S 
[Bj. At higher energies, the de Broglie wavelength becomes much less than the inter-particle 
spacing and the effects of coherent scattering are no longer important. In this limit, the bi¬ 
nary scattering approximation is recovered, i.e., E —>■ cr and pf —>■ pf. pf is more commonly 
known as the momentum transfer collision frequency, p^, which is associated with the mo¬ 
mentum transfer cross section, cr^. Similarly, Pm = kf is known as the effective momentum 
transfer collision frequency. We will adopt this convention in the following discussions. 

At higher fields, incoherent inelastic scattering effects, such as electronic excitations, 
need to be considered jSDiisg. By considering only a single inelastic channel, and assuming 
neutral particles are in the ground state, the excitation collision operator is 
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Jr Ui) = ru) 


(U + Ui), 1 = 0, 


( 10 ) 


0 , 


I > 1 . 


where Uj is the threshold energy associated with the excitation collision frequency 

III. SOLUTION TECHNIQUE 

Boltzmann’s equation is a non-linear integro-differential equation involving three spatial 
dimensions, three velocity dimensions and time. The Boltzmann equation consists of two 
parts; an advective component (in phase space) and a component representing collisions. It 
is a formidable task to solve the Boltzmann equation numerically, using a single numerical 
scheme for both components and a single time-stepping method. Because of the complexity, 
we choose to replace the task of solving the full Boltzmann equation by the task of solving 
the configuration-space transport, the energy-space transport and the contributions due to 
collisions separately, then combining the results in a manner that appropriately approximates 
the full solution. This can be achieved via the technique known as operator splitting [53] . 

To this end, the Legendre polynomial expansion of Boltzmann’s equation in plane parallel 
geometry given in equation Q can be represented as 


+ Szifi) + Suifi) — 0, 


( 11 ) 


where 



( 12 ) 



SM.) = 


p=±l 


A. Operator splitting 

The simplest method of operator splitting, and the method employed in this paper, is 
Lie-Trotter splitting [MIES], which employs two separate operators, e.g. Sz and Su, in a 
sequential order. If 


-^ + Sz{f) + Suif) = 0, 


(14) 
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then the Lie-Trotter algorithm is 


(15) 

(16) 


dt 


+ SziD = 0, with t e [r, and f* (t”) 


/ in , 


df* 

dt 


+ Su{f*) = 0, 


withie [(“,!"+'] and/*(i") 


/•(t”+‘) . 


SO that where t” and are successive times. This simple method can 

be shown to be only accurate to first order in time, and there are many other methods avail¬ 
able that offer higher order accuracy and often include additional advantageous properties 
[SSHni]. The major reason for this particular choice of operator splitting algorithm is that if 
Sz is treated in an explicit manner, and Su is treated in an implicit manner, then the result 
is essentially the Douglas class of the Alternating Direction Implicit schemes |62[ [63] , which 
is particularly successful at accurately capturing the steady-state solution. Accurately and 
consistently determining the steady-state solution can be a problem for general operator 
splitting methods [M] . 

The isolation of the configuration-space dependence to the operator Sz makes this par¬ 
ticular scheme an example of dimensional splitting. We can now investigate how we treat 
the configuration-space advection and energy-space advection and collision components nu¬ 
merically in detail. 


B. Configuration-space advection 


The operator involving the configuration-space dependence, Sz, is given by equation (12), 


which represents a coupled homogeneous advection equation. As there are no derivatives 
of U present in Sz, the configuration-space dependence can be solved independently for 
different values of U which is huge simplification when a discretization in energy space is 
used. The coupled advection equation can be simplified as follows: 




(17) 


which can be written in matrix form, 




( 18 ) 





where f = [/o, /i,and 


' 0 A<+' 

A<i“' 0 aS+’ 


( 19 ) 


^(+) 

_ 1 

^max -L 


0 


^(+) 

_ 1 

^max -L 


a\~^ 0 

I- •'max -I 

By letting A = RAR“^, where A is a matrix of eigenvalues of A on the diagonal, and R 
are the associated eigenvectors, then 


+ A^g = 0, (20) 

where g = R“^f, which now represents a set of uncoupled, homogeneous advection equations. 
It follows from the method of characteristics ra, that 


g (*. 2) = g (0, J - Af). 


( 21 ) 


Even in this extremely simple form, the solution can be troublesome. When discretized, the 
set of values z — At are unlikely to align with existing .2 values, and hence some form of 
interpolation is required. It can be shown that linear interpolation is equivalent to a hrst 
order upwind hnite volume method scheme |66]. First order methods have the advantage 
of being well behaved and can be used to conserve mass etc. with no unwanted, unphysical 
oscillations, but have the disadvantage of introducing extra numerical diffusion, particularly 
around regions of sharp variation EH. Higher order methods perform better at controlling 
unwanted diffusion but can lead to problematic, oscillatory and unphysical solutions. Rather 
than straightforward interpolation, we choose to employ a variation of a technique well 
known in fluid transport, the SHASTA algorithm of Boris and Book [68] . The SHASTA 
algorithm approach, termed flux-corrected transport (FCT), leads to a class of Eulerian 
finite-difference algorithms which strictly enforce the non-negative property of realistic mass 
and energy densities. As a result, steep gradients and shocks can be handled particularly well, 
which is a useful property when modelling transport under non-hydrodynamic conditions. 
A FCT algorithm consists conceptually of two major stages, a transport or convective stage, 
followed by an anti-diffusive or corrective stage. 

We employ a simplified version of the full FCT algorithm to numerically approximate 
g (0, 2 : — At). Let us consider the evolution of g (t, z) for a single A, i.e., g {t, z; A), over a time 
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interval of At, with a uniform configuration-space mesh with spacing Az. By discretizing in 
this way, Zj = jAz for j = 1, 2 . .., — 1, and = tn + At. The algorithm is as follows: 

1. Shift: The elements of g (t, z; A) are shifted to the node closest to Z — /3, where /? = ;^A. 
This may result in an ‘overshoot’, but we can then propagate the shifted solution 
(in step 2) either forwards or backwards in time as appropriate. The purpose of this 
step is to overcome time step limitations due to the Courant-Friedrichs-Levy (CFL) 
condition [ 62 ], which allows us to choose arbitrary time step sizes with respect to the 
configuration-space convergence (sufficiently small time steps are still necessary for the 
operator splitting accuracy etc.). By shifting to the nearest node, the CFL condition 

\d\ = I; |A| < 1 (22) 

for the remaining advection is always satisfied. 


2. Advection with additional diffusion: The advection algorithm employed is given by 


9 ^ = 


9 W+1 


9j-i) + (^7 + (fi'i+i - + g]_i ), 


(23) 


where = g {tn+i, Zj), and 
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(24) 


is the additional numerical diffusion. The dimensionless advancement P' = P — [/3] 
accounts for the shift that has been applied in step 1. Note that /5' may be opposite 
in sign to /3, which corresponds to an overshoot in step 1. However, this does not 


adversely affect the procedure. If 7 = 0, then equation (23) is the well known Lax- 
Wendroff scheme lETI, which is accurate to second order. Historically, the inclusion 
of an extra diffusion term, 7 , has been used to ensure that a density function (i.e. a 
function that is non-negative by definition) remains positive, which is unconditionally 
enforced everywhere if 7 = y. In our case, the g^ include contributions from //>i, 
which are expected to be negative in some regions of space. However, the presence of 
7 ensures the stability of which can be defined by the requirement that Ag'^~^^ < 


msix(Agj_i, Agj, Ag^_^i) where Ag^ = g'J^i — gj. When the solution g^ is sharply 
varying or, in the extreme case, a discontinuity, the additional diffusion is necessary 


to suppress unphysical oscillatory behaviour in gj 


n+1 
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3. Anti-diffusion: An ‘anti-diffusion’ step is employed to reduce the extra numerical dif¬ 
fusion introduced in (23) i.e., 



The inclusion of this extra diffusion in step 2 assures that the solution is positive and 
physically realistic, and the straightforward application of step 3 undoes this which can 
re-introduce a negative solution. Boris and Book [HH] suggested modifying the removal 
of the erroneous diffusion by just enough to maintain positivity, in a non-linear way 
(note that they worked with non-negative densities, as we have remarked on above in 
step 2). This is an early example and precursor of the modern technique of flux limiting 
[ZOMZl- In this work the full anti-diffusion step is applied in general, except in regions 
where a sharp variation or discontinuity is known a priori (e.g. configuration-space 
boundaries), in which case no anti-diffusion is applied. Unphysical oscillations can now 
occur, but we have found that for the situations considered they are negligibly small. 
The natural extension is to include flux limiting to prevent this unphysical behavior 
but this introduces extra computational complexity. The anti-diffusion step could also 
be solved implicitly rather than explicitly, but we found that this had no significant 
impact on the results. 


It should be noted that the shift step can be performed after the advection and anti-diffusion 
stages with no change in the result. We have assumed that the boundaries are absorbing, 
in that the elements oi g{t, z) that move outside the computational domain are lost, and no 
information is introduced from outside the domain. Although perfectly absorbing boundaries 
are notoriously difficult to implement numerically, in our calculations we avoid this problem 
by keeping the swarm density negligible at the simulation edges, through the use of an 


adaptive mesh, see Section HIE, In practice we pre-calculate a transformation matrix (for a 
given set of parameters) which combines the above three steps for each of the grid energies. 


C. Energy-space advection and collisions 


A major advantage of splitting the Boltzmann equation operator according to equa¬ 


tions (12)-(13) is that Su is then the familiar, spatially homogeneous Boltzmann equation. 
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There is much literature on solving this equation, and we use the approach developed previ¬ 
ously [I117SHZ7] to perform the numerical discretization and time-step. A full description of 
the numerical solution of the process is given in m, and we will briefly summarize it here. 


The equation we need to solve is equation (13). The time dimension is discretized with a first 


order implicit Euler method, which has been chosen for its good stability properties. Similar 
to the work of Winkler and collaborators [78H80] . we employ a finite difference method to 
discretize the system of ODE’s at centred points using a centred difference scheme, i.e.. 


df{U,t) 


dU 




u, 


i+l/2 


u. 


i+1 


Ui 


f{Ui+i/2,t) = 




(26) 

(27) 


so that equation (13) evaluated at i -|- 1/2 becomes, 


Su {fl)\i+l/2 = Jl (/Oli+1/2 + X] 

^ / P=±l 


^ 7 + 1/2 




u. 


i+1 


Ui 


p 


(I + + /,+p(f/,.t) 


-K 


+1/2 y 


(28) 


Although a general form can be constructed for an arbitrary grid, the simplest case is for 
evenly spaced points, i.e. 

Ui = iAU for 0 < i < n[/, (29) 

where A{7 is a constant. By discretizing at the midpoint of the two solution nodes results 
in a system of linear equations that is under-determined. The extra information is naturally 
provided by boundary conditions which are appended to the system. 

These boundary conditions have been analyzed by Winkler and collaborators [THHEO] who 
investigated the multi-term, even order approximation, and discovered that the general solu¬ 
tion of the steady-state hierarchy contains ^ (/max + 1) non-singular and | (/max + 1) singular 
fundamental solutions when U approaches infinity, and the physically relevant solution has 
to be sought within the non-singular set of fundamental solutions. The boundary conditions 
necessary for the determination of the non-singular physically relevant solution are ca 
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fi{U = 0) = 0 for odd I, 
fi{U = Uoo) = 0 for even /, 
fi{U > Uoo) = 0 for all I, 


(30) 


where Uoo represents a snfficiently large energy. In practice, Uoo has to be determined in prior 
calculation, and is chosen such that the value of foiUoo) is less than 10“^° of the maximum 
value of /o- 

D. Green’s function solution 

In our formalism and associated code, we solve for the Boltzmann equation Green’s 
function 


Cfl = <S(j - Z|,)(S(i - (o), 


(31) 


where 


^ ^ p=±l 


11 /HP+ j, (f.) 


(32) 


for / = 0,1, 2 ,..., oo. The Green’s function solution, fi, can then be used to find the solution 
of the more general space-time Boltzmann equation, i.e. 


Cfi = S{z,t), 


where S{z,t) is a source term, then 


fi {U, z,t) = / dto / dzo fi (U, z- zo,t- to) S (zq, to) • 


(33) 


(34) 


We do this by choosing an initial distribution in configuration-space that is a good approx¬ 
imation to a delta-function, which, for this study, is a narrow Gaussian, 

da{z) = —exp , (35) 

where a is a parameter controlling the width of the Gaussian, representing the temporal- 
spatial relaxation profile of a single pulse centred on Zq and released at to- In fhe limit of 
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a ^ oo, 5a{z) ^ 5{z). The formalism is quite general, enabling the treatment of various ex¬ 
periments (e.g. Pulsed Townsend (PT), Steady-State Townsend (SST) and other drift tube 
conhgurations [SI] - detailed in Section [III F ), as well as various source and spatial/energy 
space/temporal distributions, through a single solution. This approach extends the func¬ 
tionality and accuracy of Boltzmann equation solutions to those routinely achieved by Monte 
Carlo simulations musaiHai. 


E. Numerical considerations and adaptive meshing 


The matrix system of linear equations that result from the discretization of the Legendre- 
decomposed Boltzmann equation in energy- and configuration-space at each time step are of 
the size {Imax + 1)) x iP'z'nu ilmax + 1)), where and njj are the number of nodes in 

conhguration and energy space respectively. Due to the discretization schemes, the matrix is 
sparse and sparse techniques are employed to exploit this property. Each of these parameters 
are free to be increased until some convergence criterion is met. It should be noted that, 
although the two-term approximation (/max = 1) has been used extensively, it is well known 
that it can be insufhcient in many situations M- 

In order to model the spatio-temporal relaxation of a narrow Gaussian source distribution 
in conhguration-space with a distribution of energies as computationally efhciently as possi¬ 
ble, we have developed a configuration-space node-mesh that adaptively follows the size of 
the distribution throughout the simulation. In this way a small configuration-space window 
is used around the original narrow Gaussian source which can then be sufficiently resolved 
with a small n^. As the initial pulse drifts and diffuses, a small amount of information 
reaches and then leaks out of the window boundaries. Before the amount of information 
lost to the system exceeds some small tolerance, the window is extended and the solution 
at the previous time-step calculated on the new configuration-space mesh. We have found 
that the most convenient way to quantify the amount of information on the boundary is by 
the relative number density, and impose the condition that when 



n{zL or ZR,t') 
f dz n{z, t') 


> lo■^ 


(36) 


then the configuration-space window is doubled (while the number of nodes is kept the 
same). Here to is the time of the last window adjustment, zl and zr are locations of the left 
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and right configuration-space boundaries respectively. The choice to extend the window by 
doubling is to make it so that the new mesh lines up exactly with nodes of the old mesh, 
hence requiring no interpolation. The accuracy of the modified Lax-Wendroff scheme used 
to model the configuration-space advection [67] is related to the parameter /5 = ;^A, hence 
by doubling Az after a re-adjustment, the value of At can also be doubled. This effectively 
allows us to use smaller time steps when our solution is sharp and diffusing quickly, and larger 
time steps once the solution has spread out and is varying less quickly. A maximum value 
for the time step size still needs to be enforced however, since with bigger time step sizes 
less mixing between the configuration-space and energy-space components of the operator 
splitting occurs, leading to errors. 

There is one extra complication to be discussed. Since the boundaries are absorbing, 
when they are re-adjusted, the number density profiles (and distribution functions) drop 
directly from the built-up value at the previous boundaries location to zero in a single 
Az, which can lead to problematic, unphysical, oscillatory solutions when treated with the 
method described in Section |IIIB[ In order to combat this, we simply apply the procedure 
without the final anti-diffusion step for a small amount of time on the edge and in the newly 
opened regions. The extra diffusion added ensures that the solution remains positive and 
give physical results, which, after a small amount of time, ensures that the profiles decrease 
smoothly to zero at the boundary. After this short correction time, we again apply the full 
procedure. By not removing the added extra diffusion we have increased the overall diffusion, 
but since it is only applied for a small time and to a region where there is necessarily only 
a small proportion of particles, this does not significantly affect the transport profiles. 


F. Transport properties 

The cross-sections and collision operator terms represent the microscopic picture of elec¬ 
tron interactions with the medium. The macroscopic picture, e.g. transport properties that 
represent experimental measurables, are obtained as averages of certain quantities with re¬ 
spect to the distribution function, /. Among the transport properties of interest in the 
current manuscript are the number density, n, particle flux, T, and average energy, e, of the 
electron swarm, which can be calculated via 
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(37) 



(38) 


(39) 


Likewise, we can sample the traditional hydrodynamic transport coefficients in this non¬ 
hydrodynamic framework, e.g. the drift velocity, W, and the (longitudinal) diffusion coeffi¬ 
cient, Di: 



(40) 


2 



DAt) = --r — TV ( [ dz z^n 

2dt 


(41) 


where N(t) is the total number of particles: 



(42) 


When the above properties are calculated from the Green’s function solution, which corre¬ 
sponds to a simulation of a PT experiment, then the transport properties for other experi¬ 
mental systems can also be calculated in a straightforward manner. In this work we are also 
interested in the results of a SST simulation, for which there have been previous calculations 
performed for benchmark systems. Similar to [2SIEZII17], the SST transport properties can 
be determined from the Green’s function transport properties via 



(43) 


(44) 


(45) 


(46) 


In practice the upper limit of the integrals is not oo, but a sufficiently long time for the SST 
transport properties to have converged over the 2 : range considered. 
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G. Reduced variables 


Henceforth, it is convenient to work with rescaled reduced variables. In particular, the 
space and time variations will be presented as functions of 


2 ;* = noaoz, 



(47) 

(48) 


where (Tq = 10“^° m^. Likewise, the electric field dependence arises through the reduced 
electric field E/uq in units of Townsend (1 Td = 10“^^Vm^). By presenting results in this 
manner scales out the Hq dependence, and hence allows comparisons between the dilute gas 
phase and the liquid/dense gas phase, to give a true reflection of the impact of coherent and 
other scattering effects. 


IV. ELECTRON TRANSPORT IN A MODIFIED FERGUS-YEVICK HARD- 
SPHERE BENCHMARK LIQUID MODEL 

In order to investigate the effects of medium structure on charged particle transport, 
a model for the structure function is required. One such model, frequently employed in 
the literature, is that of a structure for a system of hard-sphere potentials obtained by 
applying the Percus-Yevick approximation as a closure to the Ornstein-Zernike equation, 
which yields a pair-correlation function (HU [HS]. The static structure factor is found via 
the Fourier transform of the pair-correlation function, the angle-integration of which is used 
directly in the numerical simulations. In particular, we use the model of Percus and Yevick 
with the Verlet-Weiss correction [861 IH7] to better emulate the structure of a real liquid (see 
|83] for details). The volume fraction parameter, $, specifies how tightly packed the hard 
spheres in the medium are. It can be written in terms of the hard-sphere radius r and 
the neutral number density, Hq, as d> = ^nr^rio. Smaller volume fractions indicate a larger 
interparticle spacing, and vice versa. We have modeled systems with a range of densities, 
from d> ~ 0, which approximates a dilute gas, to d> = 0.4, which states that 40% of the 
volume is excluded by hard-sphere potentials of the neutral molecules. Differences in the 
results highlight the importance of coherent elastic scattering. 

The remaining details required of the benchmark hard-sphere model implemented for 
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electron sized particles are 




^exc 


$ 


6A^ 

{ 0, U <2 eV 

0.1 A^ U>2 eV 

0, 0.2, 0.3, 0.4, 


E/tiq = 3 Td, 
mo = 4 amu. 
To = 0 K. 


(49) 


A step-like inelastic process has been included in addition to the standard Percus-Yevick 


hard sphere benchmark system in model (49). The inelastic channel introduces a peri¬ 
odic oscillatory non-hydrodynamic behaviour, similar to those observed in the well-known 
Frank-Hertz experiment, and can hence determine whether the numerical code is accurately 
capturing the non-hydrodynamic phenomena. Figure highlights the variation of the mo¬ 


mentum transfer cross-section with $, evaluated using the cross-section in (49) together 


with the static structure factor from |83]. At high energies, coherent scattering effects are 
suppressed, and the various momentum transfer cross sections converge on the dilute gas 
case (corresponding to d> = 0) . 
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Figure 1: The variation of the elastic momentum-transfer cross-section including structure 


with energy for model (49) for various volume fractions, 


The source distribution is given by 


f(U,z,0} = AfuiU}fAz), 


(SO) 


where fz{z) is a narrow Gaussian in configuration-space, i.e., 


fz{z) = 


exp 


1 / ^ 


(51) 


A^oV^ ^ ^ 2 V A^^n 

(we take Azq = 0.1), while fu{U) corresponds to drifted Maxwellian distribution with 
T = 10^ K, and W = 10^ ms“^E, i.e., 


/ m \(3/2) 


(v-W)' 


2KT 


(52) 


and A is a normalization constant such that J {U, z, 0) dU = 1. 


A. Transport coefficients in the long-time limit 


The asymptotic values of the drift velocity and diffusion coefficient calculated from the 
spatial moments (|40j) and (41) respectively using the full non-hydrodynamic code are dis¬ 


played in Table |I] for various volume fractions. Here we compare these values with those 
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determined from a purely hydrodynamic formalism and associated code [U [ZZ|. The first 
order hydrodynamic transport coefficients, i.e., the mean energy and drift velocity, agree to 
within 0 . 2 % with the asymptotic non-hydrodynamic values for the volume fractions consid¬ 
ered. The hydrodynamic and non-hydrodynamic calculations of the longitudinal diffusion 
coefficient agree to within 0.7%. As the volume fraction increases, both the mean energy, 
drift velocity and diffusion coefficient increase monotonically, a consequence of the coherent 
scattering, where at low energies, increasing volume fractions leads to decreasing structure 
factors at low Ak, and hence decreased momentum-transfer cross-sections. For further dis¬ 
cussion on the physical variation of the hydrodynamic transport coefficients with volume 
fraction the reader is referred to [501 [831 EH] • 


Table I: Comparison of the transport quantities calculated from non-hydrodynamic (first 


row) and time asymptotic hydrodynamic (second row) formalisms for model (49) at various 
volume fractions d>. 


4> 

6 

W 

noDi 


[eV] 

[10^ ms“^] 

[1024m-is-i] 

0 

0.8335 

1.385 

2.386 


0.8337 

1.385 

2.387 

0.2 

0.9765 

3.397 

6.333 


0.9772 

3.391 

6.328 

0.3 

1.080 

5.929 

11.22 


1.080 

5.921 

11.24 

0.4 

1.233 

10.52 

19.51 


1.233 

10.51 

19.63 


B. Space-time evolution of the phase-space distribution and its velocity moments 

In Figure the space-time evolution of the /o and fi velocity distribution function compo¬ 
nents are compared for $ = 0 and d* = 0.4 at three different times. The space-time evolution 
of the integral moments of %, electron density n{z,t), and velocity moment of /i, the flux 
F, are displayed in Figure 


§ 


The timescale for variation of % is governed by 


mo 


-1 
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and hence there is no explicit d> dependence in the timescale, however differences arise due to 
the implicit energy dependence in the collision frequency (which does depend on d>) and the 
coupling to higher order moments with different timescales. The timescale for variation of 
fi on the other hand is governed by 1/“^, which has an explicit $ dependence. The timescale 
for momentum exchange is significantly decreased for increasing $ at low energies, as shown 
in Figure however they approach the same value at higher energies. We will show that 
this is reflected in the evolution of the profiles. 

At small times (e.g. t* = 0.2), there are only small differences in the /o contours between 
the two volume fractions, and this is also highlighted in the density n{z,t). At higher 
energies (> 5-6 eV) there are also very little differences in the fi contours (reflecting the 
similarity in the momentum relaxation times at these energies), however at low energies, the 
d> = 0.4 contours for fi are significantly displaced in both energy and configuration-space 
relative to the $ = 0 case. This indicates significantly higher advective and diffusive fluxes 
in this energy regime at this time, which is evidenced in the flux profiles of Figure]^ Given 
the sharp initial pulse with large spatial gradients, we observe large positive and negative 
diffusive fluxes, along with a large positive advective contribution. 

At larger times, the /o and fi contours in the $ = 0.4 case depart significantly from the 
d> = 0 case, initially in the low energy regime and then finally the entire energy regime as the 
higher energy electrons relax from the initial condition. The peaks in each of the distribution 
components at larger times for $ = 0.4 case are significantly displaced in the z-direction 
from the $ = 0 case. This is reflected in both the density and flux profiles at larger times, 
which highlight the enhanced drift and diffusion due to the reduced momentum transfer 
cross-section associated with coherent scattering for this model and field. Interestingly, at 
sufficiently long times, the $ = 0.4 contours have predominantly positive values, and only 
very small negative excursions at low energies, in contrast to the $ = 0 contours. At 
these times, the flux is positive over the swarm indicating that the advective contribution 
dominates the diffusion contribution, since the density gradients are much more rapidly 
dissipated in the d> = 0.4 case, see Figure]^ 

Strikingly, both the $ = 0 and $ = 0.4 contours for both /o and fi demonstrate peri¬ 
odic structures in both configuration space and in energy space at sufficiently long times 
and sufficiently downstream from the source. The periodic structures manifest themselves 
earlier for the d> = 0.4 case. These are the well known Franck-Hertz oscillations I2ni[2ii. A 
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simplistic picture of this non-hydrodynamic phenomena is that the electrons in the swarm 
are being repeatedly accelerated by the electric field to an energy above the inelastic pro¬ 
cess threshold whereby they undergo an inelastic collision losing their energy. This simple 
physics is evidenced in the /o and fi distributions. By integrating over the energy to obtain 
the density and flux, shown in Figure much of the periodic structures observed in the dis¬ 
tribution function is masked, however some non-Gaussian spatial structure is still observed. 


We will explore the d>-dependence of the wavelengths of oscillations further in Section IV C 


C. Steady-state Townsend configuration 


The solution detailed in Section IVB|is essentially equivalent to solving for the Boltzmann 


equation Green’s function for the model (49). A strict validation of this approach and 


associated numerical code is to be able to reproduce the Steady-State Townsend (SST) 
transport properties from the Green’s function solution, as described in Section |IIID[ The 


average energy (46) and average velocity (45) for SST simulations of various volume fractions 


are shown in Figure]^ In the spatially asymptotic regime, the average energy and the average 
velocity are equal to the hydrodynamic and pulsed-Townsend values shown in Table It can 
be seen that the SST properties demonstrate damped spatially periodic structures similar to 
those observed in the Frank-Hertz experiment and other investigations [23 [2Zl EH IHMI]. 
They are a manifestation of the energy and space periodic structures in the distribution 
function components, and in the spatially periodic structures in the density and flux profiles 
of Figure By assuming the elastic scattering is weak, the width between the peaks in 
the transport property profiles. A, is directly related to the threshold energy of the inelastic 
process, Uj in eV, via 1021 

^ ' ■ (53) 


(0-^)eV/Td (-®'/'^o)Td 


where the reduced electric field is in Townsend (Td). For model (49), the theoretical spacing 
is 6.6. In Figure it is possible to see that there are variations in the wavelength of 
the spatial structures with d>, as well as significant differences in the decay rates of the 
oscillation amplitudes. For •F = 0, the wavelength is approximately 8.24 ± 0.02 and this 
decreases to 6.67 ± 0.02 for $ = 0.4. The differences arise explicitly due the differences in 
the elastic momentum transfer cross-section, as well as implicit variations associated with 
the modification to the swarm’s energy with d>. For $ = 0.4, the momentum transfer cross- 
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Figure 3: Temporal evolution of the spatially varying density and flux for model (491 with $ = 0 (dashed lines) and $ = 0.4 
(solid lines). The four columns represent the times, t* = 0.2, 2, 20, and 200 respectively. 



























section for elastic scattering is significantly reduced compared to the $ = 0 case. Hence, 
the randomizing collisions that dampen the oscillations [27] are reduced for $ = 0.4 as 
compared to other d>, and the variation of damping with $ then follows. Likewise, it should 
not be surprising that the wavelength for the $ = 0.4 case is closest to the analytic value 


of (53), since the reduced momentum transfer associated with the $ = 0.4, more closely 
approximates the weak elastic scattering assumption used in deriving it. 


We must also point out that the validity of these profiles are dependent on the discretiza¬ 
tion of the distributions in configuration-space. If the spatial discretization is of the same 
order as the Frank-Hertz wavelength, then it will be very difficult to resolve these features in 
the distributions and consequently the time-averaged profiles. Of course, our initial choice 
for the discretization is small enough to easily resolve these features, but as the simulation 
progresses and the distribution diffuses, our adaptive mesh will increase in range and also 
increase the spatial discretization step size. After a point, the coarseness of the discretiza¬ 
tion causes the distribution to slowly lose its features, which is visible in the time-averaged 
quantities by the suppression of the amplitude of the oscillations. In our simulations we 
expect our results for z* ^ 30 deviate from the true spatially dependent steady-state val¬ 
ues, however the fully relaxed values agree with the hydrodynamic values. It is simple to 
address this issue by increasing the number of points in configuration-space but this is also 
significantly more computationally intensive. 
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Figure 4: Spatial variation of the average energy and average velocity under SST conditions 
for model (49) with various volume fractions 




V. SPATIO-TEMPORAL RELAXATION OF ELECTRONS IN LIQUID ARGON 

Electron transport in liquid argon is an essential component in the function of Liquid 
Argon Time Projection Chambers (LArTPC) which are currently being used for high energy 
particle detection Ionized electrons in liquid argon originating from the high energy 
particles are accelerated under the action of an electric field to generate a current and 
consequently reconstruct the path of the high energy particle. Typically these chambers 
operate with electric field strengths of less than 500 kV/cm. The aim of this component is 
to follow the spatio-temporal evolution of these ionized electrons in liquid argon, relevant 
to the operation of these detectors. Foxe et ah [HS] have measured the energy distribution 
of the electrons ionized by high energy particles in liquid argon, and have shown that the 
majority of the ionized electrons have energies below 1 eV. Consequently in this study we 
employ an initial source energy-distribution that is constant in energy space up to 1 eV, i.e., 

fu{U) = CU-^/^eiU-ley), (54) 

where 0 is the Heaviside step function, and U is in eV and (7 is a normalisation constant. 
The mean energy of this distribution is 0.5 eV. The swarm is released from a narrow Gaussian 
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in configuration-space, 


so that the full initial phase-space distribution is f(U,z,0) = Afu(U)fz{z), where ^ is a 
normalization constant such that J {U,z,0) dU = 1. For argon, we take = 10, a 

larger initial spread than for the Percus-Yevick model, reflecting the smaller cross-sections 
of argon. 


A. Cross-sections, potentials and screening 

In a recent paper [T] we investigated the modifications required to treat transport of 
electrons in dense gaseous and liquid argon, with our simulations focused purely on the 
hydrodynamic regime. The work followed closely the approach of Lekner and Cohen ©El 
updated using modern scattering theory techniques, theories, potentials etc. The potentials, 
screening factor, and cross-sections derived in ^ are used once again here, and we briefly 
summarize the process. The gas-phase elastic scattering cross-sections were calculated by 
solving the Dirac-Fock scattering equations isi, and were shown to give good agreement 
with both beam scattering cross-sections and swarm experimental drift velocities and char¬ 
acteristic energies. As the medium’s density increased, density effects became important. 
The effects of the high density of the liquid are included in our calculations through several 
modifications of the gas-phase scattering properties, and are dependent on the liquid argon 
pair-correlation function and associated static structure factor ra. The first modification 
is to account for the screening of a single induced atomic dipole by the induced dipoles of 
all other atoms. In a high density medium, the potential produced by the induced multipole 
moments is of a sufficiently long enough range to interact with induced multipole moments 
by other atoms in the bulk. The second modification is to construct an effective potential in 
which there are contributions from both the target atom and the surrounding bulk, the latter 
of which is done as an ensemble average. The momentum transfer cross-sections calculated 
from the dilute gaseous and liquid argon potentials are shown in Figure It is significant 
to note the absence of the Ramsauer minimum in the liquid-phase cross-section. 
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Figure 5: The momentum transfer cross-sections in the gas-phase (dashed line) and liquid- 
phase (solid line) for electrons in argon [1]. 



B. Results 


To consider conditions representative of those in liquid state particle detectors, we simu¬ 
late electron transport in liquid argon under the following conditions: 


Ejn^ = 2.5 X 10“^ Td, 
T = 85 K, 
mo = 40 amu. 


(56) 


The reduced field is equivalent to 500 kV/cm with a density corresponding to liquid argon, 
no = 0.0213 A For this reduced electric field and source distribution, given in (54)-(55), 


the electron swarm energies are generally well below the first inelastic channel threshold 
energy (8.9 eV), so that there is no inelastic channel operative, and hence the periodic 
spatial structures observed in the Percus-Yevick hard-sphere liquid model above are not 
present. 

The relaxation of the /o distribution function component are compared for the gas and 
liquid phases at three different times in Figure]^ At t* = 1, there are only small differences 
between the contours reflecting similar energy relaxation rates between the two phases ini¬ 
tially. At P = 10, a bulge is beginning to develop in the gas-phase contour in the energy 
region between 0.1 — 0.5 eV, which corresponds to the Ramsauer minimum in the gas-phase 
momentum transfer cross-section. In this region, the gas-phase momentum transfer cross- 
section dips below the liquid cross-section, which has resulted in this enhancement of the 
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diffusive flux in this range. At higher energies the liquid cross-section is less than the gas- 
phase cross-section, which has resulted in enhanced diffusive flux. At t* = 100 these effects 
are even more pronounced. 

In Figure 1^ the /i component contours for the gas and liquid phases of argon are compared 
for the same three times. At the first time, t* = 1, there is already significant differences 
in the /i contours, with the largest change occurring in the Ramsauer minimum range in 
the gas-phase case. This highlights again the difference in the timescales of the energy and 
momentum relaxation between the two phases. As time increases, greater differences develop 
between the fi contours particularly around the Ramsauer minimum and at the high energy 
range for the reasons previously discussed. 

The number density as a function of time is shown in Figure The behavior of the 
number density profiles is consistent with the behaviour of the /o and /i profiles. At = 1 
there is no noticeable difference in the two number density profiles. At later times it is clear 
that, despite the Ramsauer minimum in the gas-phase, the liquid-phase experiences the 
greater diffusion rate overall. For the electric held considered and initial source distribution, 
the average drift velocity for both the gas and liquid phases is small compared to the diffusion 
rates. 





Figure 7: Temporal evolution of the number density profiles for gas-phase (dashed lines) and 
liquid-phase (solid lines) argon. The three columns represent the times, t* = 1, 10, and 100 
respectively. 
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VI. CONCLUSION 


We have developed a full multi-term, space-time dependent solution of the electron Boltz¬ 
mann equation in gases and liquids capable of modeling non-hydrodynamic conditions. The 
flexibility of the algorithm lies in solving the Boltzmann equation’s Green’s function, knowl¬ 
edge of which allows one to construct the solution for other experimental configurations 
e.g. the SST experiment and similar applications. Operator splitting has been employed 
to efficiently evolve the energy-space and configuration-space components individually with 
tailored numerical schemes. 

The theory and associated code was first applied to a simple hard-sphere benchmark 
model liquid, where structure effects were simulated by the Percus-Yevick structure factor 
as a function of the volume fraction, ff>. The inclusion of an inelastic channel was a key 
test of the algorithm’s ability to reproduce non-hydrodynamic phenomena. Periodic spatial 
structures developed in the space-time and steady-state profiles for the distribution function 
components and associated transport properties, the periodicity of which is directly related 
to the threshold energy of the inelastic process. We observed that these periodic structures 
arose on shorter times scales when coherent scattering effects became important. The steady- 
state profiles constructed for various volume fractions also reproduced the non-hydrodynamic 
oscillatory structures expected. The asymptotic transport coefficients calculated from the 
non-hydrodynamic solution of Boltzmann’s equation were also shown to be consistent with 
the values calculated from a hydrodynamic solution of Boltzmann’s equation. 

Finally, the cross-sections calculated in [1] were used to investigate the spatio-temporal 
evolution of electrons in gas-phase and liquid-phase argon. The two momentum-transfer 
cross-sections feature different qualitative and quantitative behaviours. Striking differences 
in the evolution of the components of the phase-space distribution were apparent, reflecting 
the differences in the gas-phase and liquid-phase cross-sections, particularly the absence of 
a Ramsauer minimum in the liquid-phase. This highlights the problems associated with 
treating liquid systems as gaseous systems with increased density, which has implications 
for various applications including liquid argon time projection chambers. 
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